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We present a new 3D SPH code which solves the general relativistic field+hydrodynamics equa- 
tions in the conformally flat approximation. Several test cases are considered to test different aspects 
of the code. We finally apply then the code to the coalescence of a neutron star binary system. The 
neutron stars are modeled by a polytropic equation of state (EoS) with adiabatic indices F = 2.0, 
r = 2.6 and F = 3.0. We calculate the gravitational wave signals, luminosities and frequency spec- 
tra by employing the quadrupole approximation for emission and back reaction in the slow motion 
limit. In addition, we consider the amount of ejected mass. 
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[ I. INTRODUCTION 

Coalescing binary neutron stars (NS) are expected to be one of the most powerful sources of gravitational radiation 
^ ' for detectors hke LIGO Ij], VIRGO ||, GEO600 § and TAMA @. In contrast to electromagnetic radiation, gravi- 
' tational waves (GW) are hardly perturbed during their journey from their emission site to their detector. Since the 
' GW emission is sensitive to mass distributions it provides additional information on the astrophysical object that is 
[ not available from electromagnetic radiation. To extract a signal from the noisy output of the interferometers, precise 

theoretical templates of the waveforms are needed. 
' Merging NS are also of interest to the nuclear astrophysics community since they are a promising site for the forma- 
, tion of the heaviest elements via the rapid neutron capture process (r-process). The decompression of neutron-rich 
■ material through tidal disruption has first been proposed as r-process-scenario by Lattimer |^ ^ and has been later 
— , , investigated by several groups |Q, ||, ^, |l^, [ill, Apart from that, the coalescence of two NS is at the heart of 
^ ^ several models for the still poorly understood Gamma-Ray Bursts (e.g. |l^, 13 1). 

I Since 1990, various groups have performed hydrodynamical calculations of merging NS. The first work, using Newto- 
^ ' nian gravity and a polytropic EoS was done by Oohara & Nakamura ([p^ and references therein) using an Eulerian 
Pt*] finite difference scheme. The problem has been further explored by Davies et al. |]l5|, Zhuge ||l^, |l^ and Rasio & 
Shapiro ( [p^ and references therein) using smoothed particle hydrodynamics. Improvements towards a full general 
relativistic implementation were achieved by Ayal et al. [|l9j and Faber et al. ||2^, |l] using the Post-Newtonian 
^ j approximation and Wilson et al. |p2| using the conformally flat approximation. The first fully relativistic calculations 
• have been carried out by Shibata & Uryii |2^, Hj . The microphysical aspects have been investigated by RufTert 
et al. (1^, ^ and references therein) and Rosswog et al. (|2^ and references therein), both using a realistic Nuclear 
EoS. For a review on the topic, see Q. 

The merger phase of a Neutron star binary represents only the last phase of a evolution sequence beginning with a 
long inspiral phase during which the binary evolves on nearly stationary orbits. As energy is radiated away due to 
GW, the binary separation slowly decreases down to the innermost stable circular orbit (ISGO) where hydrodynamical 
instabilities due to tidal forces set in. During the following merger phase, the two stars are disrupted and a new central 
object, either a NS or a black hole (BH) is formed in the center. This object finally settles down to equilibrium in a 
ring-down phase. Since the time scales for each phase are drastically different from each other, different computational 
approaches have to be chosen for each phase. 

The inspiral phase can be approximated to high precision by a pointmass-Post- Newtonian (FN) approximation jsj, ^ . 
This approximation is valid as long as finite-size and strong-field effects are negligible. If this is not the case, i.e. for 
very close orbits, a more exact calculation is provided by the quasi-equilibrium approximation which is based on 
quasi-equilibrium sequences in the conformally flat approximation. Very close to the ISGO this approximation breaks 
also down and dynamical calculations have to be carried out. The outcome of the merger does not only depend of the 
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GR effects but also on intrinsic properties of the neutron stars (NS) like the equation of state (EoS). The appropri- 
ate tool to treat this problem is a fully relativistic magnetohydrodynamic simulation including all the microphysics 
involved. First successes have recently been reported in implementing the GR component of the problem |5| . 

However, these calculations still remain very complicated due to the nature of the non-linearity of the field equations 
and the enormous demand in computational resources. 

Beside the PN approximation to GR there exists also the approximation using the conformally flat condition (CFC). 
Both are of about the same computational complexity. The PN approximation, which is a weak-field, slow-motion 
expansion of the field equations does not describe the problem adequately ^ in first order (IPN) since the 
expansion parameters within the NS have values of ^ 0.4 and the convergence is generally rather slow. It should be 



mentioned, however, that there has been progress in the binary-pointmass case |34 with Pade-approximants. On the 
other side, the CFC approximation can also describe strong field scenarios in special cases. The static Oppenheimer- 
Volkoff (OV) profile is reproduced exactly, uniformly rotating polytropes are described up to a precision of ~5% 
For a binary near the ISCO, estimate the magnitude of the non-diagonal spatial metric components (which are 
neglected by the CFC approximation) to ~ 0.02 for a compactness of M/R = 0.14 and ~ 0.05 for M/R = 0.19, their 
influence for the angular velocity lo and the angular momentum J to 0(10"^) for the density p, the conformal factor 
■0 and the mass M to 0(10-3). This systematical error will certainly increase in the merger phase, but the quality of 
the CFC-approximation is still not fully explored yet. 

The first hydrodynamical calculations using the CFC approximation have been carried out by Wilson with the 
irritating result of a possible collapse of both NS prior to merging. After correcting an error in the original formu- 
lation Q the method is well established today to investigate quasi-equilibrium configurations. New, independent 
hydrodynamical calculations have never been carried out since then. Doing so, such calculations could explore the 
applicability of the CFC approximation during the merger phase and eventually establish it as an economical approx- 
imation method to GR. As there exist other, much larger uncertainties on the numerical and microphysical side of 
the problem, a robust and economical simulation code based on a well-explored approximative method could be a 
step ahead. 

With this new code, using the CFC approximation, we therefore follow two goals: On one side, a fast and stable code 
to simulate compact objects is made available, on the other side, we want to explore the applicability of the CFC 
approximation during the merger phase by comparing our simulations to full GR results. 

The paper is organized as follows. We present in section II the numerical method with the adaption of the smoothed 



particle hydrodynamics (SPH) method to the CFC approximation. In section III we test and calibrate our code by 
applying it to two test scenarios that have previously been investigated with different methods. The central part of 
the paper, the simulation of the dynamical coalescence of two neutron stars is presented in section Finally, we 
conclude the paper in section 



II. NUMERICAL METHOD 



We use a Lagrangian particle method, the smoothed particle hydrodynamics (SPH), to implement the hydrodynamic 
part of the system of equations (A5- A8) while the field equations ( A10| - A17) are evaluated with a multigrid algorithm. 



The system of equations is summarized in appendix For further references to SPH, see ^ 



A. Basic SPH equations 



Since the coordinate conserved density p* in (A4) obeys a Newtonian-like continuity equation (A5) we define p* in 
SPH by 

pl = Y^m,Wab (2.1) 

b 

where mf, is the rest mass of particle h and Wab = W{\fa — rj], h) denotes the weight given by the standard spherical 
spline Kernel function W{r, h). 

The pressure gradient in the momentum equation ( |A6| ) is calculated in a manner analogous to standard SPH, replacing 
P by p*, 

1 v4 f Pb , Pa 
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The other terms in the momentum equatio n dep end on values of the gravitational potentials and their derivatives. 
Their evaluation is de scri bed in the section ( 11 C ). 

The energy equation (A8) consists of Newtonian- like part which can be evaluated in the same manner as in standard 
SPH 



dt 



IE 



nib 



Pa 



Pb 



PaPl PbPb 



Vb)VWab 



^-^ln(az.V^)a. 
Pa at 



(2.3) 



(2.4) 



The total time derivative of \n{au'^ip^)a in the second term, which arises due to the temporal change of the metric, 
is evaluated by second order finite differencing in time. The hydrodynamic evolution is governed by the equation of 
state 



whereas initial data is generated by 
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(2.5) 
(2.6) 



where k is the polytropic constant and F the adiabatic index. 

Since physical units enter only through the constant k, the whole problem can be written in a non-dimensional way 
by the following transformation [p9| 



(2.7) 



Therefore, k can be removed completely from the equations. Once results are obtained they can be rescaled to physical 
quantities via (2. 7). To integrate the system in time, a second order Runge-Kutta scheme with adaptive timestep is 
used. The timestep control mechanism is chosen as in m with Cat = 0.8. Trial runs with smaller timesteps lead 
to similar results. The simulations are carried out on a standard 500 MHz Athlon workstation and last about three 
weeks, where the evaluation of the field equation in the nested grid setup takes most of the time. 
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B. Artificial Viscosity 

We also include an artificial viscosity (AV) scheme to treat shocks. The standard scheme proposed in |Q is known 
to give good results in shock situations but to introduce spurious viscosity forces in shear flows. We therefore choose 
a scheme [|3| which assigns each particle a its own time dependent viscosity coefficient aa[t). The coefficient is 
determined by integrating the simple differential equation 

d eta — amin „ /o o\ 

-aa = \- Sa (2.8) 



dt 



Ta 



where Sa = max(0, —{Vv)a)- The first term of the RHS causes aa to decay on a timescale Ta towards a minimum 

amin- The second term leads to a rapid increase of only in case of a shoc k ( i.e. (Vw)a largely negative). 

To incorporate the AV scheme in a relativistic consistent manner we follow |44| by adding the viscous pressure term 

Qa 

(Vw)a + /3a^a(Vw)^)a (Vu)a > , , 

^'^ ~ [0 else ^ ■ ' 

to the physical fluid pressure, where Pa = '^cta, Wa — 1 +Pa/Pa + fa is the enthalpy, ha the smoothing length and Ca 
is the speed of sound. 

We directly use the SPH estimate for the velocity divergence 

{^V)a « - 5Z "^biVa ' Vb)^ aWab (2.10) 



C. Evaluating the field equations 



All six elliptic field equations for aijj, and x (AlC, Al^ ,Al6 A17) are solved on an overlaid gravity grid via a 

multigrid solver (e.g. |45[). 

The transfer of all the quantities involved - hydrodynamical quantities as well as gravitational fields - between SPH 
particles and the grid points is done by assigning the SPH-interpolated values to the grid points and by interpolating 
the grid values to the SPH particles using a TSC scheme known from particle mesh codes 

To take into account the matter density distribution during the late stages of a neutron star merger (rapidly rotating 
central object plus low density material surrounding it), we additionally implement a nested grid algorithm which 
allows us to enlarge the computational domain while keeping the resolution in a selected subdomain. To extend the 
multigrid algorithm to the nested grid setup, we used the Multi Level Adaptive Technique (MLAT) (see e.g. |^ ). 
Beyond these two grids, the potentials are extrapolated to the few particles which flow over the coarse grid. We use 
a quadrupole expansion around both individual stars for ip 
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where 



S^{x)d^x 
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denote the monopole and the quadrupole of the potentials' source. The index n runs over both stars since we expand 
around each star separately. A similar expression holds for aip. These expansions follow from the field eqns. (A13) 
and dAlOD . 

For the shift-vector potentials, we make use of the imposed boundary conditions (A21) - (A24) 
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(2.12) 
(2.13) 
(2.14) 
(2.15) 



where is the particle position of particle a and its position projected onto the grid boundary. 



D. Backreaction Force and Gravitational Wave extraction 



To implement the radiation reaction and the GW extraction formulae (A25,A32) we consider the slow- motion 
quadrupole moment in SPH approximation 
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To calculate gravitational wave sig nal and luminosity, we proceed as follows. The second term in the brackets in 
much smaller than the first one |60|| and therefore neglected. To obtain we take numerical derivatives of the above 
expression. This quantity is stored as a function of time and then smoothed and derivated numerically again in a 
postprocessing step to obtain the third time derivatives. With this data we are able to calculate the gravitational 
wave luminosity and the angular momentum loss without any further numerical approximation after the simulation. 
In addition, we determine the frequency spectrum using 



^(/) = |(W)/^(|/.+ (/)P 



where and (^x) denote the angle-averaged waveforms |If 



14 , 
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(2.19) 
(2.20) 
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and and {h^) are their Fourier transforms. 

A simpler approach needs to be chosen for the fifth quadrupole derivative in the radiation reaction force. We assume 
the Quadrupole to be constant in the binary's corotation frame. Therefore the only time dependence arises from the 
rotation of the coordinate frames which is easy to evaluate and differentiate analytically. This leads to fairly good 
results as long as we have a binary with two distinct companion stars. However, for a typical merger event, this 
approximation has shown to underestimate the angular momentum loss by ~ 40% during GW luminosity peak (see 
Fig. Blanchet et al. proposed a radiation reaction scheme involving only third order derivatives. However, 
the quadrupole derivatives and the radiation reaction force formula are only approximate to Newtonian order which 
rules out at least a direct implementation of this scheme. 
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- dj/dt using apprdximated derivatives 
dJ/dt using exact derivatives (postprccessing) 




FIG. 1: Comparison of tlie angular momentum loss with approximated derivatives and with exact numerical derivatives by 
postprocessing of the quadrupoles. Shown for model A (see section IV) 



III. CODE CALIBRATION 



In this section, we calibrate our code by examining three important test cases for which data from different ap- 
proaches is available. The first one is a single rotating polytrope, the second test case is a corotating quasi-equilibrium 
binary system and the last one tests the gravitational wave luminosity of binary systems. 
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A. Rapidly Rotating Polytrope 



Axisymmetric configurations in general relativity have been treated in full GR as well as with a CFC approximation 
scheme [^5[ We compare the results obtained against those from |Q. To achieve this, we prepare a F = 3, Mq = 
0.1484 - neutron star with an density profile corresponding to a static spherical star and an initial velocity field 

v = (ixf. (3.1) 

The star is then relaxed into an uniformly rotating equilibrium state with given angular velocity Q using the relaxation 
scheme described in appendix The values for fl have been chosen to fit the values in table I of . 
We measure the stars' angular momentum, the gravitational mass, the central energy density Cc — p(l + e)|c and the 
ratio of rotational to gravitational energy T/W. 



TABLE I: For each value of the angular velocity fl, we list in the upper line the central energy density Cc, the gravitational 
mass Mq, the angular momentum and the ratio of rotational over gravitational energy T/W . In the lower line, we list the 
values from [p5| 



n 




Mc 


J. 


T/W 


0.0000 


0.960 


0.1233 


0.0 


0.0 


0.0000 


1.000 


0.1232 


0.0 


0.0 


0.3315 


0.910 


0.1236 


0.00353 


0.0149 


0.3332 


0.903 


0.1237 


0.00356 


0.0149 


0.4872 


0.824 


0.1245 


0.00556 


0.0355 


0.4881 


0.815 


0.1245 


0.00559 


0.0358 


0.5920 


0.751 


0.1254 


0.00726 


0.0590 


0.5927 


0.736 


0.1254 


0.00730 


0.0597 


0.6630 


0.674 


0.1265 


0.00891 


0.0855 


0.6632 


0.665 


0.1265 


0.00891 


0.0858 


0.7061 


0.610 


0.1275 


0.01043 


0.1122 


0.7064 


0.600 


0.1276 


0.01049 


0.1132 


0.7250 


0.558 


0.1285 


0.01178 


0.1358 


0.7256 


0.542 


0.1287 


0.01208 


0.1413 





We find generally an agreement within 0.1% in Mq, ^ 1% in Jz and ~ 1% in T/W. The more sensitive local 
quantity Cc agrees to 3%. 

In the last example we faced problems to relax the configuration because the chosen angular velocity is almost at the 
mass-shedding limit. 

We used a model with 9939 particles and a (65x65x65) grid. The grid spacing has been chosen to cover the star by 
~ 20 mesh cells. 



B. Corotating Quasi-Equilibrium Binary systems 



As a second test case, we have chosen to investigate the properties of corotating quasi-equilibrium binary systems, 
a system which lies already quite close to our final merger application. We compare our results to ]49[| . 
To construct the system, we set up a pair of Oppenheimer Volkoff stars of equal mass on a given orbit with a fixed 
binary separation and a guess angular velocity. The relaxation scheme described in appendix ^ drives the system 
towards a corotating quasi-equilibrium binary system. We store angular velocity, angular momentum, gravitational 
mass and the relative separation z — Vin/rout, i-e the ratio of the innermost to the outermost point of the stars. 
Four relaxation processes with varying binary separation have been performed. The results are listed in table || 
together with results from Q . 

We find a relative discrepancy of ~ 0.1% for Mq, ^ 0.5% for Jz and ^ 4% for z. Again, z as a local quantity is more 
sensitive. The binary has been modeled with 2x25173 particles and we used a (65x65x65) grid such that one star 
has been covered by ^ 18 grid cells. 
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TABLE II; For each value oi Q, we list in the upper line the relative separation z, half of the gravitational mass Mq and half 
of the angular momentum . In the lower line, we list interpolated values from 





z 


Ma 


J. 






n 14000 


041 R 


0.116 


0.202 


0.14086 


0.0418 


0.123 


0.173 


0.14090 


0.0418 


0.123 


0.169 


0.14085 


0.0417 


0.130 


0.135 


0.14090 


0.0419 


0.130 


0.135 


0.14084 


0.0416 


0.134 


0.120 


0.14092 


0.0419 


0.134 


0.115 


0.14084 


0.0416 



C. Gravitational Wave Luminosity from Quasi-Equilibrium systems 



The radiation reaction and GW extraction mechanism in our code is independent from the CFG im plementat ion 
and needs to be calibrated separately. To do so we compare our slow- motion quadrupole formalism ( A28| , A31 ) to 
the results of | |5^ which relies basically on a quasi-equilibrium approach in the GFG approximation but with a more 
accurate GW extraction scheme at the grid boundary. 

First, we construct an equilibrium configuration using the same method as above. Then, the radiation reaction scheme 
is switched on and the GW-luminosity is measured. Using the same binary system that was used in |]5^ , we calculated 
four values in a range which was accessible both the QE approach and our code. The results are plotted in Fig. 0. 
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FIG. 2: Gravitational wave luminosity of a quasi-equilibrium system with A/o=0.1 and r=2. Duez et al. extract the GW 
signal at the boundary of the computational grid while in this work the signal is extracted using the slow-motion quadrupole 
formula. 



Mo reover, we test internal consistency by comparing the angular m omen tum loss due to backreaction force 
( A25| ) and angular momentum loss predicted by t he q uadrupole formula ( A3C ) . Both curves use the approximated 
quadrupole derivatives. Fig. ^ shows that ansatz (A25) is able to reproduce the slow- motion quadrupole formalism. 



IV. BINARY NEUTRON STAR MERGER 



Having done the tests, we apply the code to the merger of two NS in a binary system which are modeled as 
simple polytropes. Three models are considered with a different adiabatic index each time. We analyze morphology, 
gravitational wave signal, gravitational luminosity spectrum and amount of ejected material depending on the stiffness 
of the EoS. 
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dj/dt due to backreaction force 
dj/dt predicted by quadrupole formula 




FIG. 3: Internal consistency test of the backre actio n scheme. Compared are the angular momentum loss due to backreaction 
and the prediction of the quadrupole formula ( A3C ) . Both calculations use the approximated quadrupole derivatives. 



A. Initial Setup 



For a model with a given adiabatic index F, we have basically two free parameters to choose, the rest mass Mq and, 
if we want to work in physical quantities, the polytropic constant k. We assign to the models with a stiff EoS a larger 
rest mass since as a stiff EoS allows a much larger maximal mass than a soft one. The polytropic constant k is chosen 
to match the radii of the three models to i?o — 9.6km. To choose a optimal initial binary separation, the ISCO has 
been determined experimentally by relaxing the binary systems in the mutual field. By varying the binary separation 
is steps of 0.5km, we find the closest binary which is still stable against tidal disruption but whose separation is close 
to the ISCO to save computational resources. 

For all three models, we use 50346 SPH particles in total for the hydrodynamics and a (65 x 65 x 65) grid to solve 
the field equations. As the central high density region forms, we switch to a nested grid configuration (see Fig. Q). 
The relevant parameters are summarized in table IIL 



TABLE III: Input parameters for models A - C. Listed are the adiabatic index F, polytropic constant k, the ratio of binary 
separation at the start of simulation to the radius of a spherical star diniuai/ Rq, the rest mass of one star, half of the 
total gravitational mass, the ratio of rest mass of one star to the maximal mass of the OV star of the corresponding model 
Cmass = Mo/M™"^ 



Model 


F 




dinitial/ Ro 


Mo 


Mg/2 


Cmass 


A 


2.0 


7.032 X 10' 


3.02 


1.34 


1.235 


0.89 


B 


2.6 


5.600 X 10^ 


3.13 


1.60 


1.421 


0.70 


C 


3.0 


1.071 X 10" 


3.47 


1.74 


1.555 


0.65 



In the following discussion we measure distance in units of Rq, density in units of po, the central density of the 
OV profile, time in units of P, the rotation period of the binary where d = 3i?o, mass and energy in units of Mq and 
the GW- frequency is measured in units of /o = 2/P. Moreover, we choose the origin of our time axis at the GW 
luminosity peak. To allow the reader a conversion into physical units, we summarize the values of the normalization 
constants in table |^ 



TABLE IV: Normalization constants used in the discussion 



Model 


-Ro [km] 


Po [gem ^] 


P [ms] 


Mo [Mq] 


fo [Hz] 


A 


9.6 


1.30x10''' 


2.051 


1.34 


975.27 


B 


9.6 


8.66x10" 


1.955 


1.60 


1023.27 


C 


9.6 


7.64x10" 


1.935 


1.74 


1033.58 
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FIG. 4: Gravity grid setup during the simulation. The top figure sketches the single-grid setup with a (65 x 65 x 65) grid during 
the pre-merger phase. The bottom figure shows the two-level-nested grid configuration with two (65 x 65 x 65) grids. 

B. Morphology 

The simulation is started from the quasi-equilibrium configuration outside the ISCO. Since this setup contains no 
inward radial velocity, the binary needs initiaUy about half a revolution to start with the orbital decay. The stars 
then slowly approach each other and merge within about three orbital revolutions. 

Shortly after the two stellar cores begin to merge, the gravitational wave luminosity reaches its maximum, followed by 
the minimum as the two cores begin to build up the central object. At the same time, the system starts to shed mass 
through the outer Lagrangian points and builds up spiral arms around the newly forming central object. Contrary 
to Newtonian simulations |2^, the spiral arms do not extend far out into space, but get closely wrapped around 
and finally sink back into the disk of matter around the merger site. In the post-merger phase, all models lead to a 
differentially rotating central object which first bounces and then slowly contracts (see Fig. The deviation from 
axisymmetry disappears within 3 (model A) to 7 (model C) revolutions. 




FIG. 5: Normalized density distribution in the orbital plane of Model A. The contour lines range from 10 ^ to lO" with a 
logarithmic spacing of 0.25 dex. 
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C. Evolution of Angular Momentum and Gravitational Mass 



An important test of every hydrodynamical code is the check for conserved quantities. In our case, we store angular 
momentum and gravitational mass which would be conserved to IPN order in the absence of backreaction forces 
In Fig. o we plot the two quantities versus time for the three mo dels. Likewise, we plot the same quantities. 



but with the integrated angular momentum and energy loss (A3C,A31) subtracted to follow the net conservation 
without backreaction. 

Concerning angular momentum, we see in all models a loss of 8% arising from the backreaction force. In the 
post-merger stage, this force is artificially suppressed in our calculations due to our approximation of the quadrupole 
derivatives (see Fig. ^) and the loss in the angular momentum of ~ 0.5% is mainly due to numerical errors caused 
by the limited resolution of the gravity grid. The corrected angular momenta are conserved within 1% during the 
whole merger process. This is obviously well within the expected IPN accuracy range of {M/Kf' ^ 5xlG~^. 
A similar behaviour is observed for the gravitational masses. After a decrease of ~ 0.4% due to GW radiation, the 
masses of all models decrease by ~ 0.2% due to the same numerical errors. The small jumps are due to the resizing 
of the gravity grid during the run as the expression for Mq contains the non-compact term KijK^K The corrected 
masses vary within ^ 0.2 — 0.3% over the whole process. 




r^2.o 
r=2.6 
r=3.o 

r=2.0, corrected 
r=2.6, corrected 
r-3.0, corrected 



FIG. 6: Plotted is the temporal evolution of the total angular momentum and gravitational mass. The angular momentum 
is normalized to the initial value, the gravitational mass to total rest mass. Also indicated are the same corrected quantities 
but with the integrated angular momentum and energy loss subtracted to show the net conservation. The angular momentum 
decrease due to gravitational radiation is ~ 8% but the corrected angular momenta are conserved within ~ 1%. The jumps in 
the gravitational mass occur due to resizing of the gravity grid setup. 



D. Maximal densities and Black hole formation 



For corotating quasi-equilibrium sequences, a decrease of the maximal density with decreasing binary distance is 
predicted We therefore spend special care to the maximal density evolution during inspiral and relax the 

binary as carefully as possible prior to the simulation. In Fig. 0, we plot pmax versus time. During the inspiral phase, 
we see for all models a nearly constant evolution with a small increase of 1% for model C, followed by the sharp 
decline at t/P ~ as the stars get disrupted. This effect is comparable to the numerical accuracy of the code. 
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At the end of the merger phase, a dense central object builds up. It contracts but soon bounces back due to pressure 
and centrifugal forces. This bounce can be seen in all three models by a peak in central density, i.e. none of the 
models does immediately collapse to a BH. Afterwards, an oscillating central object with rising central density forms. 
However the grid resolution to calculate the gravity during the ringdown stage is not sufficient. Test calculations 
with enhanced resolution have shown increased self-gravity forces and higher central densities. We therefore consider 
the central densities in the present simulation as a lower limit. 

1.6 



1.4 



1.2 

1 

0.8 
0.6 

-2-10123 
t/P 

FIG. 7: Time evolution of the maximal density in the three models. The values are normalized to the central densities po of 
the OV star corresponding to the model. The soft model starts with the lowest maximal density of the three since the tidal 
distortion is largest in this case. The maximum drops to a minimum as the cores begin to merge. Also visible in all models is 
a central bounce just after the merger phase. In the ringdown phase of the central object, the maximal density steadily rises 
without convergence, towards a possible later collapse. 

We can also consider the dimensionless rotation parameter a = J^/Mq of the system. If this parameter is less that 
unity, in most cases, except for very stiff EoSs, a black hole will be formed. This parameter drops below unity for all 
three models after the merger phase. 

Model A can be compared to model (C2) of |25| . The parameters ainiuai, the rotation parameter at the beginning 
of the simulation, and Cmass, the ratio of the rest mass of one star to the maximal mass for a. spherical star allowed 
by the EoS agree within 2% and we therefore expect a qualitative agreement. However, |^ see a collapse to a 
BH on a dynamical timescale without a bounce contrary to this work. We see three possible reasons for the different 
outcome. 

• The numerical resolution is comparable in both models, but not before we switch to the nested grid mode. This 
happens just before the bounce, at about t/P ~ O.I. Therefore, it is well possible, that we underestimate the 
central densities in the moment when an eventual collapse builds up. 

• The slightly different dynamical evolution during the merger phase when the CFC approximation deviates most 
from full GR could have a direct influence. 

• Since differential rotation plays an essential role in the stability of supramassive stars, we investigated the 
rotational pattern of the central object. We fit the angular velocities lUi to a typical aj(r) = ti;c/(I + ''ck;^"^) " 
law 1^ to get a quantitative estimate and obtain a value for A which is a measure for the amount of differential 
rotation. For model A, the result is A ~ I.7i?o, i-e A ~ A/Re ~ 1.1 if we take Re — 1.5Ro for the equatorial 
radius of the central object. But an object with a value of A~^ ~ 0.9 could well support a rest mass larger 
than 3Mq[|51||. This is a strong argument against the collapse in our calculations, a difference in the numerical 
viscosity might therefore explain the different outcome. 



E. Ejected Matter 

With respect to NS mergers being a promising site for r-process, a crucial question is the amount of matter that 
becomes unbound and gets ejected into space. In Newtonian simulations, a small fraction at the tip of the spiral arms 



r=2.6 
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FIG. 8: Plotted are the angular velocities uii vs. the cylindrical radius at t/P = 1.19 for model A. The angular velocity is 
normalized to the angular velocity where d = 3Ro 



has enough total energy to leave the central gravitational potential. To determine the amount of unbound matter in 
the relativistic case, we need a criterion similar to the Newtonian energy-criterion. To do so, we assume the following 
conditions for a spiral arm particle: 

• Pressure forces are negligible with respect to gravitational forces (i.e. the particles move on geodesies) 

• Particles move in a stationary metric (which is produced mainly by the axisymmetric, rapidly rotating central 
object) 

Then, we can find a conserved expression CstaUonary which is similar to the Newtonian Ckin + Sint + epot-Expression 
for total energy per unit mass 

^stationary — V Ui pr -\- pr 1 (4.1) 



The derivation can be sketched as follows. Using the momentum eqn. (AC) we get 

p*^(.^u~ ) = p*^Mu.) - l(uO)p*-^^ ~ 2ap*uOa (4.2) 
at at at (|y0^2^4 

where we used dta = 0, 9t/3* = 0, dttp = 0, p = and the time derivative of the identity of = — + ^"^'^„ . On the 
other hand with the aid of the energy eqn. (A8) we get 



Adding up the two identities and using eqn. ( |A7| ) we arrive at 



d , , ™ 1 

which simplifies to 



,.^,„.„,_^.„, + 2-i, + JI^ + 2-l) = (4.4) 



To estimate how much mass will be ejected, we plot ^stationary as a function of the distance to the central object at the 
end of the simulation. There are only several distinct particles having positive energy corresponding to a cumulative 
rest mass of approximative 1.5xlO~'*M0 (9 particles), 5xlO^''^M0 (2) and 2xlO~'*M0 (5) for models A, B and C, 
respectively. 

These values can only be taken as an order of estimate since the resolution in the spiral arm region is far too low. 
Compared to model C in |l^ (Afej ^I.SxIO^^Mq) which is the Newtonian analogon to model B of this work, the 
relativistic effect of a stronger gravitational attraction clearly becomes apparent. 

It has been stated, that given the uncertainties in the merger rates, even an amount as small as a few times lO^^M© 
per event in r-process material may be an important contribution to the enrichment of the Galaxy with heavy elements 
(see Fig. 26, in ]To| , with an upper limit to the merger rate of 10~^ per year and galaxy js^). Therefore, further, 
high-resolution calculations are needed for precise answers. 
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FIG. 9: Plotted are the values of tstauonary versus the cylindrical radius r^yi = ^ -\- ]p- at tjP = 1.19 for model A. The 
particles with values below zero are considered to be bound. Only very few particles are unbound corresponding to some 
1O~*M0. The situation is similar for model B and C. 

F. Gravitational Radiation 

The gravitational radiation waveform as well as the gravitational wave frequency spectrum are expected to contain 
important information about the dynamics of the merger and the internal structure of the NS, especially the EoS. In 
the following, we compare the radiation waveform and frequency spectrum depending on the stiffness of the EoS. 
The gravitational radiation waveforms are shown in Fig. |l^. To allow for comparison, we matched the simulated 
waveform onto a pointmass result which is based b asical ly on Newtonian gravity p4|, but with the inspiral time t 



fitted to the code signal and with the quadrupoles (A28) used. Apart from a small accumulating phaseshift which 



results from the Newtonian approximation, we notice that the code signals follow remarkable closely the pointmass 
signal until only a few oscillations before merging. Partly a reason for this effect is the fact that we work in the slow 
motion limit. Only the F = 2 signal shows an exceptional behaviour as it exceeds the pointmass signal in amplitude 
and frequency at about tjP ~ —0.4, two oscillations before merging. This is caused by a dramatical acceleration of 
the orbital decay and a speedup of the binary in angular velocity. We do not see the effect in such a drastic form in 
the stiffer models B and C. The reason for the special behaviour of model A is a numerical one: the error in angular 
momentum is up to one third of the backreaction loss during the last orbit which causes the observed speedup. A 
comparative run to model A with a slightly different numerical setup shows a different angular momentum violation 
behaviour. This results in a different acceleration of the orbital decay and the angular velocity and in slightly different 
wave signals. 

A clear difference between the three models, which has also been pointed out in the Newtonian case is the 
intensity and length of the ringdown signal. The GW signal of the soft-model A shows a distinct peak before decaying 
rapidly within ~10 oscillations, whereas the stiffer models are characterized by larger peaks and a decay on a longer 
timescale which exceeds the simulation timescale in model C. Contrary to other investigations [|l^ a third peak 
cannot be seen. 

All these phenomena can also be observed in the luminosity spectrum Fig. |ri|. We indicate the GW-frequencies at 
dynamical instability and at maximal gravitational wave luminosity, jdyn and jaw-, respectively. The quadrupole 
oscillation frequency of final central object is denoted by fQpoie- In the frequency range below fdyn, the three models 
follow the pointmass behaviour dE/df{f) ~ f~^^^- Beyond this frequency, model A begins to deviate, due to the 
same numerical reason mentioned above: while the angular velocity increases, the merging process is accelerated and 
less energy (which is the time integral of the GW-luminosity) will be accumulated in one frequency interval. The 
other two models follow the pointmass template up to about few- The spectrum of the comparative run to model 
A also deviates from the pointmass spectrum at /^yn but lies above that of model A. At log(///o) — 0.45, the curve 
then bends down to the minimum. 

Since all our simulations underestimate the backreaction force in the late merger phase, the real merging process 
might exhibit a larger acceleration than the one in our simulations. Therefore, we expect that the real spectra will 
lie below ours and therefore also below the pointmass spectra in the frequency range between fdyn and few- Such a 
behaviour has been observed in Newtonian simulations [|7| . 

The ringdown signals translate into strong maxima, peaked around the quadrupole oscillation frequency fgpoie- Here 
we observe clear differences in strength between the three models. While the stiff-EoS model C peaks well above 
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the pointmass line the peak of model A is much weaker. This has already been shown in Newtonian jlj] and Post- 
Newtonian investigations. Moreover the normalized quadrupole frequency fqpoiel fo in model A is slightly higher 
than in the other two models since in model A, the soft EoS produces a much more centrally condensed central object. 
Therefore, we are able to distinguish the stiff-EoS spectrum from the soft-EoS spectrum by means of the ringdown 
part. However, as long as the numerical inaccuracies during the merger phase are still of the order of the backreaction 
terms, the exact shape of the part below few cannot be reliably determined. 

Contrary to other investigations in PN approximation we do no observe a dip between fdyn and few- The 
difference may be explained in the different pointmass waveform fitting procedure. While [s^ adds a purely Newtonian 
waveform, we modify the inspiral time r and the quadrupoles. This increases the frequency of the pointmass waveform 
so the dip is filled up. 
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FIG. 10: Gravitational waveforms of models A to C normalized to the rest mass. The time axis is normalized to the orbital 
period P at a binary distance d = SRo- The dashed lines correspond to fitted pointmass waveforms. 
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FIG. 11: Frequency spectra for the three different models. The dashed line shows the pointmass spectrum dE/df ~ f^^^^ ■ The 
frequency axis is normalized to the GW frequency /o at a binary distance d = "iRo- Model A deviates early from the pointmass 
behaviour because of numerical errors. Due to a build-up of a non-axisymmetric central object, the stifF-EoS model C produces 

a much higher quadrupolc oscillation peak than the soft-EoS model A, whose central object quickly relaxes to axisymmetry. 
Furthermore, as model A produces a more centrally condensed object, the quadrupole oscillation peak is shifted to slightly 
higher frequencies. 
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FIG. 12: Gravitational wave luminosities for the three models considered. All models show a second peak but no further peaks 
can be detected. The strength increases with higher stiffness. 



V. CONCLUSION 

We have developed and tested a 3D smoothed particle hydrodynamic code based on the conformally flat approxi- 
mation to GR with an additional approximative backreaction scheme. The gravity field equations are solved on an 
overlaid grid with a multigrid solver. The code has been applied to the neutron star merger problem using a polytropic 
equation of state. Three models have been considered, with adiabatic indices F = 2.0, F = 2.6 and F = 3.0. We start 
with quasi-equilibrium models to avoid oscillations. In all models, we find about 8% of angular momentum being 
radiated away due to backreaction. 

None of the central objects does immediately collapse to a black hole. Instead, we always find a bounce followed by 
a slow contraction. We also find a highly differentially rotating object in all three models. 

The amount of unbound mass is significantly less than obtained in Newtonian simulations. We obtain values up 
to some 10~'*Mq. This mass corresponds to only a small number of particles, therefore final conclusions cannot be 
drawn. 

We are able to distinguish the gravitational wave signal and spectrum of a stiff-EoS model from a soft-EoS model by 
means of the ringdown part because the stiff-EoS model produces a much stronger signal during this phase. In general, 
however, we are not yet able to determine the exact shape of the spectrum due to numerical errors caused by the 
limited resolution of the overlaid field solving grid and a too simple evaluation of the quadrupole derivatives around 
the gravitational radiation peak. In future work we will increase the resolution by increasing the particle number and 
using a finer field solving grid to make the results more reliable. The code shall also be generalized to nonsynchronized 
binaries, an initial condition which is likely the case for real systems. Moreover, we plan to incorporate a realistic 
EoS (e.g. 0). 
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APPENDIX A: CONFORMALLY FLAT FORMALISM 

We shortly summarize in this section the formalism on which our code is based. It follows closely the derivation in 
psf and [|9j, but is slightly adapted to a Lagrangian form suitable for SPH. Concerning indices, we use the summation 
convention where Latin indices run from 1 to 3 and Greek indices from to 3. We also set G = c = 1. 
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1. Conformally flat approximation 

The spatial part 7ij of the metric is approximated by the conformally flat form 

lij = (Ai) 

so that the ADM line element reads 

ds^ = {-a^ + PiP')dt'^ + 2Pidx'dt + i^^Sijdx'dx^ . (A2) 
Here, a is the lapse function, pi the shift vector and ijj the conformal factor. 

2. Hydrodynamic equations 

The fluid is described with a perfect fluid stress-energy tensor 

T^. = ,9«;wV+p5''" (A3) 

where p refers to the rest mass density, w = l+ p/p + eto the specific enthalpy, p to the pressure, e to the specific 
internal energy, and u'^ to the four velocity of the fluid. 
With the coordinate conserved density 

p* = pau°il)^ (A4) 

the continuity equation reads 

j^p* + p*dy = 0, (A5) 

where = u'^/vP is the physical velocity. The relativistic momentum equation reads 

4:Ui = — %oi'4}^dip - avPdia + Ujdi(3^ + '^^''^'^ diip, (A6) 
at p* ^5y0 

where Up = wUf^ is the specific momentum which connects to the physical velocity by = — /3' + ^"^p . Moreover, we 
get from the normalization condition U/^u^ = — 1 

iauy = l + ^. (A7) 
An energy conservation equation enables us to evolve the internal energy 



and finally, an equation of state 



closes the system. 



^, = _Pdy-Plln{au'^^% (A8) 
dt p pat 

P = P{p, e) (A9) 
3. Field equations 



Equations for the metric components a, /3* and V can be found from the maximal slicing assumption, the 
Hamiltonian constraint and the momentum constraint. 

The maximal slicing assumption tiKij = 0, where Kij is the extrinsic curvature, leads to a Poisson-like equation 

A(a^) = 27Ta^^{E + 2TijSijtp-'') + lai^^KijK'^ (AlO) 
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where 



and 



E — pw{au^)'^ — p 



The Hamiltonian constraint serves us an equation for the conformal factor 

8 

And, finaUy, with the momentum constraint we obtain an equation for the shift vector 



where Si = p*Uiip ^. Using the definition 



( A14) sphts into two simpler parts 



AS' = a In. ,^ 



and 



Ax = djB^. 

Furthermore, we extract the foUowing physical observables 

Mg = -2 J S^d^x, 



Mo 



(All) 
(A12) 

(A13) 

(A14) 

(A15) 

(A16) 

(A17) 

(A18) 
(A19) 
(A20) 



which correspond to the total rest mass, total angular momentum and total gravitational mass. We have to be aware, 
that these quantities are conserved only to IPN order, the order of the CFC method. Therefore, measurements of 
these quantities are only accurate to order (M/R). 



4. Boundary conditions 



To solve the equations for aip ( A10| ) and -0 (A13), we expand the source around both stars in multipoles and 
approximate the values on the outer boundary of the computational grid up to the hexadecupole order. Note that 
due to the finite grid size, it is not possible to include the contribution from KijK^^ which extends beyond the grid 
boundary. However, by comparing the contribution from KijK^^ on different grid sizes, and using its r~^-fall-off- 



dependence, which follows from the asymptotic form of the potentials and (A12), we assume the exterior contribution 
to be < 0.01 Mq. 



The same argument, on the other hand, does not hold for the shift-vector equations (A16) and (A17) since the non- 



compact-support source term dj ln(a'0 ^) {djfi^ + di(3^ — ^Sijdif]^) is not small compared to the matter source term 
and we therefore cannot calculate exactly the integrals over the source function Sf. A possible solution is to impose 
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only fall-off boundary conditions |49j (note: use a different coordinate system) which makes them independent of the 
absolute values of the exact integrals. 

S°d^x = b'= (A21) 

= (A22) 

(A23) 
(A24) 



5. Backreaction Force and Gravitational Wave Extraction 

Like the Newtonian and the PN approximation, the CFC approximation does not include gravitational radiation, 
an additional scheme to mimic the loss of energy and angular momentum and an extraction scheme of gravitational 
waves needs to be implemented. The radiation reaction force is modeled via an ansatz which generalizes the Burke- 
Thorne radiation reaction formula. The gravitational wave extraction scheme is approximated to quadrupole order. 
For the radiation reaction force, we propose the following ansatz 

where 

F'-^'^^ = -\x,Xj{Q,ji^'^ (A26) 
o 

and <7 is the active gravitational mass density in the CFC approximation 

a = T0° + T" = p(l + ^ + e){uy{l + v^) + - ^ + S^'^). (A27) 

The quadrupole finally is approximated [ p2| 

= STF |-2y S^{x)x^xj(f-x\^ . (A28) 



where S^{x) is the source term in equation (A13). We recover the quadrupole in the so-called slow-motion approx- 



imation which is valid for slowly moving strong- field sources |59 . The symbol 'STF' refers to 'symmetric trace- free 



part' which can be calculated in the case of a two- index quantity as 

STFA.j = -I- ]^Aj, - ^5^jAkk. (A29) 

This approximation is based on the idea to expand a metric perturbed by radiation into mass and current multipoles 
of the radiating source |59[| . Using the above quadrupole, we get up to quadrupole order 

= -^e»,fe(Qj™)^')(Qfe™)'''. (A30) 

Law = MG = -^(Qy)'''(Qy)^'^- (A31) 



To justify the scheme, we compare ansatz ( A25D with (A30) in Fig. || to ensure internal consistency and (A31) with 



results from in Fig. |^. The evaluation of the fifth quadrupole derivative is described in section II 
The gravitational wave signal is extracted in the TT-Gauge 

hfi^it)^lPvkiQki{t-r) (A32) 

where Pijki — (Sik — nink){6ji — rijUi) — l/2(6ij — ninj){5ki — nkUi) is the transverse-traceless projection tensor and 
n, = Xil\yi\. 

Due to parity, the mass-current quadrupole and the octupole moment vanish in the systems considered in this paper. 
Otherwise we would also have to consider these terms. 
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APPENDIX B: RELAXATION SCHEME 



In order to start the simulation with weU-defined initial conditions, we need a scheme that produces a binary 
configuration with a preselected velocity field. In our case, we considered synchronized binaries, i.e. binaries with a 
velocity field 

v = nxf (Bi) 

where |f2| is chosen such that the attracting gravity force is balanced by the centrifugal force. 
To achieve this goal we add a relaxation acceleration that drives the fluid into equilibrium 

/ = —{v-nxr) (B2) 

where is given by the balancing condition 



Here ri2 is the binary separation and Vr — v ■ Cr- We approximate v by straight forward finite differencing in time 

u v{t) - v{t - dt) 



dt 



(B4) 



Eqn. (B2) is equivalent of adding a braking force Jbr = —{v — Q x r)/Treiax- Since not but rather Ui are the 
fundamental variables, we have to update them by means of 

Ui = Ui{v) ^ {v' + f3')wu°i;'^ (B5) 

to take to changes the desired effect. 

In order to keep the binary separation at a desired value ri2, we adjust the stars' position from time to time and set 
the velocity field w to a synchronized one. Then, we also need to readjust the whole field and specific momentum 
distribution while keeping the physical velocity field constant. This is done by iteratively solving the field equations 
([A1^)-(|aT7|) and 

U^,new = (1 - 5)u^,old + Siv' + f3')wU°^P^ (B6) 

where S is typically chosen to be 0.3 to 0.5 to avoid overshooting. 
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